Numerical integration of orbital motion

Example1: Orbit around Earth (central force)

Here is a simple example of how to use the force functions (plus a few other ones) defined in Celestlab to integrate the motion.

Only the central force is considered. The trajectory is compared with a purely Keplerian trajectory.

// -----------------------------------------------------
// Utilities
// -----------------------------------------------------
// State evolution: dy/dt = f(t, y). 
// y = [pos; vel] (state vector)
// t in seconds
// model = central force
function [ydot]=fct(t, y)
  pos = y(1:3,:); 
  vel = y(4:6,:); 
  ydot = zeros(y);
  ydot(1:3,:) = vel;
  ydot(4:6,:) = CL_fo_centralAcc(pos); 
endfunction

// -----------------------------------------------------
// Initializations, integration
// -----------------------------------------------------
// initial state
t0 = 0;  
pos0 = [7000.e3; 0; 0]; 
vel0 = [0; 5.e3; 5.e3];
y0 = [pos0; vel0]; 

t = (0:180:10000); // time instants from t0 (seconds)

// integration
rtol = 1.e-12 * [1;1;1;1;1;1];
atol = 1.e-6 * [1;1;1;1.e-3;1.e-3;1.e-3]; 

y = ode(y0, t0, t, rtol, atol, fct);
pos = y(1:3,:); 
vel = y(4:6,:); 

// Comparison with Keplerian motion 
kep0 = CL_oe_car2kep(pos0, vel0); 
kep = CL_ex_kepler(0, kep0, (t-t0)/86400); 
[pos_k, vel_k] = CL_oe_kep2car(kep); 

// -----------------------------------------------------
// Plots
// -----------------------------------------------------

scf();
plot(t/86400, CL_norm(pos-pos_k),"k"); 
xtitle("Error on position (m)", "Time (days)"); 
CL_g_stdaxes()

scf();
plot(t/86400, CL_norm(vel-vel_k),"k"); 
xtitle("Error on velocity (m/s)", "Time (days)"); 
CL_g_stdaxes()

Example2: Orbit around Earth (central force + J2..J6)

Here is an example of how to use the force functions (plus a few other ones) defined in Celestlab to integrate the motion.

Initial position and velocity are initialized from mean orbital elements (Eckstein-Hechler model). The trajectory is compared with that obtained by Eckstein-Hechler. The Energy is also computed, which gives an idea of the integration accuracy.

// -----------------------------------------------------
// Utilities
// -----------------------------------------------------
// State evolution: dy/dt = f(t, y). 
// y = [pos; vel] (state vector)
// t in seconds
// model = central force + zonal terms (J2:J6)
function [ydot]=fct(t, y)
  pos = y(1:3,:); 
  vel = y(4:6,:); 
  ydot = zeros(y);
  ydot(1:3,:) = vel;
  ydot(4:6,:) = CL_fo_centralAcc(pos) + CL_fo_zonHarmAcc(pos, 2:6); 
endfunction

// Potential (for integration checking)
function [pot]=potential(pos) 
  pot = CL_fo_centralPot(pos) + CL_fo_zonHarmPot(pos, 2:6); 
endfunction

// -----------------------------------------------------
// Initializations, integration
// -----------------------------------------------------
// orbital elements = sma, ex, ey, inc, raan, psom 
t0 = 0; 
param0_mean = [7000.e3; 0; 0.001; 98*%pi/180; 0; 0]; 
param0_osc = CL_ex_propagate("eckhech", "cir", t0, param0_mean, t0, "o"); 
[pos0, vel0] = CL_oe_cir2car(param0_osc); 
y0 = [pos0; vel0]; // state vector
t = (0:180:86400); // time instants from t0 (seconds)

// integration
rtol = 1.e-12 * [1;1;1;1;1;1];
atol = 1.e-6 * [1;1;1;1.e-3;1.e-3;1.e-3]; 

y = ode(y0, t0, t, rtol, atol, fct);
pos = y(1:3,:); 
vel = y(4:6,:); 

// Conversion to mean orbital elements 
param_osc = CL_oe_car2cir(pos, vel); 
param_mean = CL_ex_meanElem("eckhech", "cir", param_osc); 

// propagation with Eckstein-Hechler 
[param_mean_eh, param_osc_eh] = CL_ex_propagate("eckhech", "cir", 0, ..
                                param0_mean, (t-t0)/86400, "mo"); 
[pos_eh, vel_eh] = CL_oe_cir2car(param_osc_eh); 

// Energy
mu = CL_dataGet("mu"); 
energy = CL_dot(vel)/2 - potential(pos);

// -----------------------------------------------------
// Plots
// -----------------------------------------------------

scf();
plot(t/86400, CL_norm(pos-pos_eh),"k"); 
xtitle("Position difference with Eckstein-Hechler J2..J6 (m)", "Time (days)"); 
CL_g_stdaxes()

scf();
plot(param_mean(2,:), param_mean(3,:), "b"); 
plot(param_mean_eh(2,:), param_mean_eh(3,:), "or"); 
xtitle("Mean eccentricity vector"); 
CL_g_stdaxes()

scf();
plot(t/86400, -mu./(2*energy) + mu./(2*energy(1)),"k"); 
xtitle("Energy variation (m)", "Time (days)"); 
CL_g_stdaxes()

Report an issue
<< Frame conversions	Cookbook	Rotations >>